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Abstract. - Incommensurate structures can be described by the Frenkel Kontorova model. 
Aubry has shown that, at a critical value K c of the coupling of the harmonic chain to an 



incommensurate periodic potential, the system displays the analyticity breaking transition 
between a sliding and pinned state. The ground state equations coincide with the standard 

■ map in non-linear dynamics, with smooth or chaotic orbits below and above K c respectively. 

For the standard map, Greene and MacKay have calculated the value ^=-971635. Conversely, 
evaluations based on the analyticity breaking of the modulation function have been performed 
for high commensurate approximants. Here we show how the modulation function of the infinite 
system can be calculated without using approximants but by Taylor expansions of increasing 
order. This approach leads to a value _K^=. 97978, implying the existence of a golden invariant 

' circle up to K' c > K c . 



> ' 

Introduction. The Aubry transition in the Frenkel Kontorova model is considered 
equivalent to the breakup of tori in the standard map Q. Based on Greene's hypothesis 
that dissolution of invariant tori can be associated with the sudden change from stability to 
instability of nearly closed orbits, an accurate evaluation of the critical coupling K c has been 
derived Besides, by a rigorous computer proof |I| an upper limit K c < 63/64 has been 

established above which no golden invariant circles can exist. Here we give a comparably 
accurate estimate of K c based on the analyticity breaking of the modulation function. Our 
approach does not make use of finite commensurate approximants, but of Taylor expansions 
up to very high order of the modulation function of infinite systems. The fact that the two 
values of K c do not coincide seems to imply violation of Greene's assumption. 
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The model. - The Frenkel Kontorova model (FKM) consists of a harmonic chain of atoms 
interacting with a potential V of period incommensurate to the average lattice spacing I. The 
total potential energy reads: 



- 1 

S = ^-(x l -a; l _i-0 2 + V(xi) (1) 



=i 



with V{x) = 77-^(1 - cos(2ttx)) (2) 
( 27I T 

where Xi is the coordinate of particle i. The usual choice is to set the period of the potential 
to unity and the average spacing to the golden mean I = r g = (\/5 — l)/2. 

The ground state positions xf can be described by the modulation function g: xf — 
il + g{il) where g{x) is periodic with the period 1 of the modulation potential. 

Correspondence with the standard map. - At equilibrium, the force on each particle has 
to vanish. The relation §^ = 2xi — x i _ 1 — x i+ i + V'(xi) — yields the standard map. 

f \ _ j f Xi \ _ ( 2x.i + V'(xi) - 

V x i ) \ Xi ~ 1 / \ Xi 

Starting from a point (xi,Xq), this relation defines iteratively a sequence of points on a torus 
(xi+i, Xi) mod I?. For K < K c subsequent points lie on smooth orbits whereas for K > K c 
all orbits become chaotic The theoretical explanation of this transition follows from the 

Kolmogorov, Arnold and Moser (KAM) theorem 0. 

The standard map and the Aubry transition can be related by considering the trajectories 
with an infinitesimal change in starting point and their relation to the sliding phonon mode 
in FKM. A small change of starting point (Sxi, Sxq) will lead to successive changes from the 
original trajectory of eq. (0), which up to first order in <5a; are: 




Sxi+i \ ( 2 + V"(xi) -1 \ / 5x t 
Sxi { 1 Oil Sn-i 



(4) 



The time dependence of phonon eigenvectors e of the FKM, ej = Xi(t) — xf ~ e lwt , satisfy: 

Cj+i = 2ei - 6i-i + (V"(xi) - u?)ti (5) 

When u> = 0, eq. (||) for 6i becomes equivalent to eq. (Q) for Sxi. If the standard map 
trajectories are unstable with respect to initial conditions, the divergence of Sx would imply 
localization of the sliding mode if lo = belongs to the spectrum Q. Since for u> = 0, 
6j = 1 + g'(il) one can easily imagine that an exponentially localized eigenvector yields a 
discontinuous modulation function g. This in turn leads to pinning of the FKM and to the 
disappearance of the sliding mode. Interestingly we have noted that in the neighborhood 
of K c a tendency to localization of the lowest frequency mode occurs, showing up in a drop 
of the participation ratio. We have conjectured that localization of the zero-frequency mode 
is the precursor of the structural transition which leads to the opening of the phonon gap. 

The critical value K c has been determined to a high precision as K c = 0.971635406 by 
MacKay by renormalization of invariant circles within Greene's hypothesis In this paper, 
we propose a route to the calculation of K c which does not need this hypothesis. 
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Modulation junctions at the Aubry transition. Aubry has discussed the analyticity 
breaking transition also by applying a perturbative approach to the equilibrium equation for 
the modulation function. 

2g(x)-g(x + l)+g(x-l) = -V(x + g(x)) (6) 

The KAM theorem could be applied to prove that there exists a continuous function, which 
obeys eq. (||), when the potential is smooth enough and K is small. We follow the empirical 
approach of ref. |9| , which gives insight in the small denominator problem. By writing V' in 
Fourier series as 

+00 

V'{x) = K V m e 2mmx (7) 



m— — 00 



and approximating V {x + g(x)) ~ V {x) in eq. (Q) we find the modulation function to be: 

g{x) = -KY^^ mx (8) 

with uj m — 2 1 sin(7r?7?i)|. Notice that this approach applies to a general periodic potential 
V. The denominators uj^ can become infinitesimally small, especially for Fibonacci numbers 
F n (F =Ft = l,Fi + + Fj_ 2 ). With the relation F n _i - F n r g = {~T g ) n+1 the divergence 
of the denominators are: 

iim 4- = A(4r +1 ( 9 ) 

For sufficiently small values of K, this divergence is limited. Aubry H gave, with a simple 
derivation, an upper bound for K above which an analytical modulation function can no 
longer exist. This value for potential (§) is An, which is much larger than the actual K c . Our 
goal is to find the value of K for which divergence of the denominators really takes place 
and compare this value with the one obtained by MacKay H . The equivalence of these two 
approaches relies on the validity of Greene's hypothesis. 

Approximations of the modulation function of infinite systems and evaluation of K c . 
We rewrite eq. (^|) with potential (||) by applying the discrete Fourier transform: 

+00 1 
g(x)= x ke 2 ™ kx with inverse: X k = / dx g(x) e - 27rlkx (10) 



yielding: 



e e y J — e e J 



(ii) 



By expanding exp(27rig(a;)) as in Q eq. ( |Tl| ) becomes: 



2 Y Ki (' 2?r ) 



in 



E 



Afej Xk 2 ■ ■ ■ Xk m 5k 1 +k 2 +...+k m +i,k 



m=0 kik2-..k,_ 

— (—l) m Xk 1 Xk a ■ ■ ■ Xk m 6k 1 +k 2 +...+k m -i,k 



(12) 
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To transform eq. (|i2|) into an expansion in K we define: 



X k = KXi 



K 2 Xl + K 3 Xl 



(13) 



yielding: 



h,k - S-i, 



4-7T m! 

m— 1 



■ E E 

nin!...iim k 1 k 2 ...k„ 



9m-t-n2 + ---+ro m + l,n 



fci fc2 



' X k™$k 1 +k 2 + ...+k m + l,k - {-^) mX kl X kl 



' X k™^k 1 +k 2 + ...+k m -l,k 



(14) 



Eq. |TJ could have been obtained directly by applying to both sides of eq. (0) the operator 



D — n\ dh 



n — Jo dx 



-2-xikx 



approximation g(x) = — ^ Sln ^ 7r ^ 4- 0(K 2 ). Since higher harmonics Xk scale with K k [|10[ 
= for |fc| > n. A manageable iterative algorithm to calculate the coefficients X^ from 
those of lower order can be derived by defining the matrix P(n, k, to) as: 



|jf=o- From eq. ( |l4] ) it is straightforward to obtain first order 



P(n, k, to) = — — p- ^ ^ ^fe^ 1 ^ 2 • ■ ■ ^™<5„ 1+ „ 2+ ... +nm: „4i+fc2+...+fc m ,fc (15) 

nin 2 ...n, T1 k 1 k 2 ...k m 

P(n, k, to) = for |fc| > rt and for m > n. Eventually we are only interested in the elements 
with to = 1 which give the Taylor- Fourier coefficients of the modulation function 

XI = (16) 

However elements of P with to > 1 are used during the calculation, allowing an effective 
iterative procedure. To this purpose we rewrite eq. (|l^) for m > 1 as: 

P(niU) ^ryr™N!_! y y 

v ' ' ' to k ™ (m-l)\ ^ 

n m km n 1 n 2 ...n m - 1 k\k 2 ...k m -\ 

X ki X k 2 ■ ■ ■ X k ™ ^5 ni+n2 + ... +rim _ 1 ^ nm 6k 1 +k 2 + ... + k m - 1 ,k-k m ^ (17) 

This leads to our final recursive relations: 

1 n— 1 

P(n,k,l) = --Lj- 2 y\p{n~l,k~l,m)~-{~l) m P{n-l,k + l,m) (18) 

TCI — 1 

^ ra — rn+1 min{n',fc+n— ri' } 

P(n,k,m) = — y y P(n',k',l)P(n- n',k- k',m- 1), (n > to > 1) 

n' — l k' — max{ — n' n+n'} 

The second relation is nothing but eq. ( |i~4| ) written in term of P for m = 1; The first one 
starts the iteration for n = 1. The third relation uses eq. ( |l7| ) to calculate elements with 
to > 1. All elements of P(n, fe, to) are real and the set of equations (|l8|) is very efficient for 
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numerical calculations. Memory is the bottleneck at very high Taylor order, since the number 
of nonzero elements up to order n, which have to be stored, increases as ~ n 3 . CPU time also 
increases with n 3 but remains moderate, a 150-th order calculation taking about 15 minutes 
on a Pentium. 

The iterative calculation of eq. ( |l8| ) has to be performed only once to get, by use of eq. ([l6|) , 
( |l3|) and (|Io|), all modulation functions of order n for K S [0 : K c ]. For K > K c we expect 
this method to break up either because eq. (|l3| ) is no longer converging for some k or because 
infinite high-frequency oscillations appear in the discrete Fourier sum. 

In fig. [I] we plot the modulation functions for different K, together with the forces f(x) = 
2g(x) — g(x — I) — g(x + I) + V'(x + g(x)), which should vanish for the exact modulation 
function. The appearance of high frequency oscillations at K = 1 > K c is clearly visible, 
yielding non- vanishing forces. In fig. || we show the Fourier coefficients up to Taylor order 
n = 150 for several values of K across the transition. Their values agree with numerical 
evaluations on commensurate approximants jy]. The behaviour as K k expected for K — > 
is still a good representation even at K = 0.5 ~ K c /2. This behaviour is found for all k and 
is not limited to k < 20 by numerical accuracy as found in . 

Inspection of the modulation function as given in fig. |l| is not enough to pinpoint K c . 
However the behaviour with increasing n of the coefficients Xj? can lead to a precise evaluation 
based on the following reasoning. By assuming that the Taylor-Fourier coefficients grow with 
a power law: 

X£~P(n,fc,l)~ (A*) n 

X: ~ P(»,n,l) ~ (An)" (19) 
both K n X™ and the sum (|l^) remain convergent up to KX < 1 whence we can estimate K c 
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Fig. 3 - \P(n,k,l)\ as function of n for different k for the FKM. The plotted curves are 
max(|P(n, k, 1)|) (the k ma x value of k corresponding to this maximum is given in the lower panel), 
\P(n, n, l)|,|P(n, 1, 1)| and \P(n, 2, 1)|. Note that when n is a Fibonacci number, the maximum is 
always at k = n. The straight line is the fit y = a\ n through P(89, 89, 1) and P(144, 144, 1), giving 
A = 1.012 => K c (Fn) ~ 0.988. Dotted lines in the lower panel correspond to the Fibonacci numbers. 



as: 



K ° = 77 — TT ( 2 °) 

In fig. || we show the elements \P(n,k, 1)| as function of n for several values of k. We 
also show max(|P(n, k, 1)|), i.e. the maximum of \P(n, k, 1)| for each n. The fc-value corre- 
sponding to this maximum, indicated as k max , is shown in the lower panel. In particular for 
n corresponding to Fibonacci numbers (34, 55 ,89, 144), k max is equal to n. Moreover, at 
these values of n we see sudden jumps due to the fact that uj^ 2 diverges like eq. (^|). Besides, 
a power law dependence of the maxima at Fibonacci numbers begins to develop for large n 
as indicated by the straight line. 

From these numerical calculations we can conclude that A„ is the dominant exponent. 
Therefore for the calculation of K c we need to calculate eq. [l8| only for k = n where it takes 
a simpler form 

^ n— 1 

P(n,n,l) = — lu- 2 V P(n- l,n- l,m) (21) 

m— 1 
_^ n— m+l 

P(n,n,m) = — N P(n' ,n' ,l)P(n — n' ,n — n' ,m — 1), (n > m > 1) (22) 

m — j 



m 

n' = l 



For this case only ^n(n + 1) elements have to be stored, and calculations up to the order 
of a few thousands can easily be achieved in order to get a good estimate of A„. In fig. [| a) 
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The line aX n corresponds to the best fit through n — 
>onds to 0.979778542, which is our best approximation 
3 as a function of F n . Note that the order is up to 
the upper limit 63/64 and MacKay's value ||. For 
and quadruple precision are included up to Fig. Only 
but not systematic. 



P{n, n, 1) is calculated up to n = F20 and found to have a power-law behaviour. As we assume 
that P(F n ,F n , 1) ~ (tz~) F ": we can define successive approximations of K c as: 



K c (F n ) 



P(F n ,F n ,l) 



P(F n -l, -Fn-lj 1) 



(23) 



In fig. |b) we show K c (F n ) up to order F 20 = 109 46, with K C (F 2Q ) = 0.979778542, which 
is close but significantly different than MacKay's value K c — 0.971635406 fl. The value of 
K c = 0.979 < 63/64 determined by our procedure seems well converged and is unlikely to 
reach MacKay's value for higher F n . We have also checked that quadruple precision does not 
change our numerical results. The discrepancy between our value of K c , which is directly 
related to non analyticity of the modulation function, and that of the transition to instability 
in the standard map questions the validity of Greene's hypothesis. 
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